# Table 5: Robustness Tests


# 1. Load Packages ----

library(lfe)
library(stargazer)
library(miceadds)
library(margins)
library(tidyverse)
library(xtable)
library(haven)

# 2. Read in Data ----

load(file = "df_voxit_individual_props.RData")

load(file = "df_voxit_individual_analysis.RData")

# 3. Prepare Data ----
voxit_individual <- voxit_individual %>%
  mutate(sum_binary1 = ifelse(utility_optc_sum>0, 1, 0),
         max_binary1 = ifelse(utility_optc_max>0, 1, 0))

# 4. Conduct Regressions ----
## 4.a Turnout ----
### i) Logit Regression (Column (1)) ----
logitx_vcov <- glm.cluster(formula = turnout/100 ~ utility_optc_sum + utility_optc_max +
                             factor(year) + factor(canton_id) + prop_voxit +
                             male + married + age + uni + knowledge + leftright +
                             initiative + initiative_counter + referendum_fak + 
                             referendum_comp ,
                           cluster = as.numeric(voxit_individual[voxit_individual$prop_voxit>1,]$datum_merge) +
                             voxit_individual[voxit_individual$prop_voxit>1,]$canton_id,
                           family = binomial(link = "logit"), 
                           data = voxit_individual[voxit_individual$prop_voxit>1,])
logitx <- glm(formula = turnout/100 ~ utility_optc_sum + utility_optc_max +
                factor(year) + factor(canton_id) + prop_voxit +
                male + married + age + uni + knowledge + leftright +
                initiative + initiative_counter + referendum_fak + 
                referendum_comp ,
              family = binomial(link = "logit"), 
              data = voxit_individual[voxit_individual$prop_voxit>1,])
margx <- summary(margins(model = logitx, vcov = logitx_vcov$vcov, 
                         variables = c("utility_optc_sum", "utility_optc_max")))

# Regression Table
# Utility sum
# p-value
margx[[2,5]] #***
U_sum <- c(paste0(round(margx[[2,2]]*100, 3),"***"))
U_sum_se <- c(round(margx[[2,3]]*100,3))
# Utility max
# p-value
margx[[1,5]] #***
U_max <- c(paste0(round(margx[[1,2]]*100, 3), "***"))
U_max_se <- c(round(margx[[1,3]]*100, 3))
Observations <- length(logitx$data$turnout)
c1 <- rbind(U_sum, U_sum_se, U_max, U_max_se, Observations)
# Save Column (1)
print(xtable(c1, digits = 3), file = "Table5_column1.tex")


### ii) Binary Utilities (Column (2)) ----

lm.bx <- felm(turnout ~ sum_binary1 + max_binary1 + prop_voxit + 
                male + married + age + uni + knowledge + leftright +
                initiative + initiative_counter + referendum_fak + referendum_comp | 
                canton_id + year | 0 | canton_id + datum_merge, 
              data = voxit_individual[voxit_individual$prop_voxit>1,])
summary(lm.bx)


# Regression table
stargazer(lm.bx, 
          type = "latex", star.cutoffs = c(0.1, 0.05, 0.01),
          star.char = c("*", "**", "***"), summary=T,
          keep = c("sum_binary1", "max_binary1"),
          covariate.labels = c("1\\{U^{sum}>0\\}", "1\\{U^{max}>0\\}"),
          df = F, dep.var.caption="Dependent variable: Turnout",
          dep.var.labels.include=F, float = F, 
          omit.table.layout ="n",
          out = "Table5_column2.tex")



## 4.b Selective Abstention ----
### i) Logit Regression (Column (3)) ----
l_e_vcov <- glm.cluster(formula = empty_ballot/100 ~ utility_optc + prop_nr_day +
                          male + married + age + uni + knowledge + leftright +
                          haven::as_factor(typex) * factor(datum_merge) + canton_id,
                        cluster = voxit_propositions[voxit_propositions$turnout==1,]$id +
                          as.numeric(voxit_propositions[voxit_propositions$turnout==1,]$datum_merge),
                        family = binomial(link = "logit"), 
                        data = voxit_propositions[voxit_propositions$turnout==1,])
l_e <- glm(formula = empty_ballot/100 ~ utility_optc + prop_nr_day +
             male + married + age + uni + knowledge + leftright +
             haven::as_factor(typex) * factor(datum_merge) + canton_id ,
           family = binomial(link = "logit"), 
           data = voxit_propositions[voxit_propositions$turnout==1,])
# CAUTION: time-consuming
margx_empty <- summary(margins(model = l_e, vcov = l_e_vcov$vcov, 
                               variables = c("utility_optc", "prop_nr_day")))


# Regression Table
# Net Benefit
# p-value
margx_empty[[2,5]] #***
Net_benefit <- c(paste0(round(margx_empty[[2,2]]*100, 3),"***"))
Net_benefit_se <- c(round(margx_empty[[2,3]]*100,3))
# Utility max
# p-value
margx_empty[[1,5]] #***
Ballot_position <- c(paste0(round(margx_empty[[1,2]]*100, 3), "***"))
Ballot_position_se <- c(round(margx_empty[[1,3]]*100, 3))
Observations <- length(l_e$data$empty_ballot)
c3 <- rbind(Net_benefit, Net_benefit_se, Ballot_position, Ballot_position_se, Observations)
# Save Column (3)
print(xtable(c3, digits = 3), file = "Table5_column3.tex")



### ii) Excluding Observations (Column (4)) ----
# Excluding observations that make the identifying assumption implausible
lm.eb_rob <- felm(empty_ballot ~ utility_optc + prop_nr_day +
                    male + married + age + uni + knowledge + leftright | 
                    haven::as_factor(typex) * factor(datum_merge) + canton_id | 0 | id + datum_merge, 
                  data = voxit_propositions[voxit_propositions$turnout==1& voxit_propositions$problem_date!=1,])
summary(lm.eb_rob)

# Regression Table Column (4)
stargazer(lm.eb_rob,
          keep = c("utility_optc", "prop_nr_day"),
          covariate.labels = c("Net benefit", "Ballot position"),
          out = "Table5_column4.tex")



